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Abstract. We study the order-disorder transition of the two dimensional interacting 
monomer-dimer model (IMD) which has two symmetric absorbing states. To be self- 
contained, we first estimate numerically the dynamic exponent z of the two dimensional 
Ising model. From the relaxation dynamics of the magnetization at the critical point, 
we obtain fil{vz) = 0.057 650(12), or z = 2.168 26(45), where /3 = | and u = 1 
are exactly known exponents. We, then, compare the critical relaxation of the order 
parameter at the transition point of the IMD with that of the Ising model. We found 
that the critical relaxation exponent /3/{vz) is in good agreement with the Ising model, 
unlike the recent claim by Nam et al [JSTAT ( 2014 ), P0801I]. We also claim that 
the Binder cumulant is not an efficient quantity to locate the order-disorder transition 
point of the model with two symmetric absorbing states. 
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In two dimensions, systems with two symmetric absorbing states generally exhibit two 
phase transitions m [21 El SI El E], One is the symmetry breaking order-disorder transition 
(SBODT) and the other is the absorbing phase transition (APT). When these two 
transitions coincide [7], the critical relaxation dynamics shows logarithmic behavior 
like the voter model [8] and, in this context, the universality class to which models 
with two symmetric absorbing states belong is termed as the generalized voter (GV) 
class [siiniiio]. 

By capturing main feature of models with two symmetric absorbing states, A1 
Hammal et al |2] suggested a representative Langevin equation of the GV class, which 
resembles the model A (according to the classification scheme of [11]) with a single 
component. Because the critical behavior of the model A with a single component 
order parameter is robust against various nonequilibrium perturbations [121 ESI ITTj|l. 
the SBODT, once occurring at a distinct point from the APT point, is expected to 
share criticality with the Ising model. This indeed was confirmed numerically for the 
two-dimensional interacting monomers (2DIM) model [6j. 

Recent numerical study of the two dimensional interacting monomer-dimer model 
(IMD) which also has two symmetric absorbing states, however, challenged this general 
picture and the SBODT of the IMD was claimed not to be the same as the Ising-type 
phase transition [16]. Nam et al [16] argued that non-Ising criticality can be originated 
from interfacial fluctuation of the third state, so-called ‘dimer’ state which is absent in 
the 2DIM. Hence more extensive numerical study for the IMD is desired to settle down 
the issue of the universality class. The purpose of this paper is to clarify the universality 
class of the IMD by extensive Monte Garlo simulations. 

This paper is organized as follows: In section [2l the dynamic rules of the two 
dimensional IMD are explained. Section El presents simulation results. Since the 
relaxation dynamics of the Ising model at the critical point will play an important role 
in this paper, we present simulation results about the dynamic exponent of the Ising 
model in section 13.11 Then simulation results for the IMD are presented in section 13.21 
Section 0] summarizes and concludes this work. 

2. Two dimensional Interacting Monomer-Dimer Model 

As a variant of a catalytic surface reaction model proposed in m, the IMD was first 
introduced as a one dimensional model with two symmetric absorbing states [18] and 
the generalization to two dimensions was introduced and studied in [5] [16] . This section 
explains the dynamics of the two dimensional IMD through a simulation algorithm and 
introduces quantities we are interested in. Let us hrst denote the index of each lattice 
point by n = (i, j) (i, j = 1,2,..., L). Periodic boundary conditions are assumed. For 

f For the model A with multi-component order parameter, certain nonequilibrium perturbations are 
relevant m- 
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convenience, a site (i, j) will be called an even (odd) site if i + j is even (odd). 

Each site is one of three states; A-occupied, i?-occupied, and vacant. Each site n 
is assigned a random variable which takes one of three values 1, 0, and —1 if the 
site is occupied by A (cr^j = 1), occupied by B = 0), or vacant (cr^j = — 1). By 
[Bf^ is denoted the number of A’s [BA] in the nearest neighbors of site n. If an = —1 
(vacant) and Afi ^ 4, the site n is referred to as an active site. If there is no active site, 
no change of conhgurations is allowed. In this sense, a conhguration without any active 
site is absorbing and there are two absorbing states; all even sites are occupied by A 
and all odd sites are vacant, and vice versaj^ 

Now we are ready for explaining an algorithm for simulations. Let us assume that 
there are active sites at time t. We choose one of active sites at random with equal 
probability. Assume that site n = (i,j) is selected. With probability p a ‘monomer’ A 
attempts to adsorb at this site (monomer-event) and with probability 1 — p a ‘dimer’ 
BB attempts to adsorb at site n and one of nearest neighbors of n (dimer-event). For 
a dimer-event, the nearest neighbor site fh can be either [i -|- l,j) or {i,j -|- 1) with 
probability Note that even if we choose m among four nearest neighbors with equal 
probability, the result is the same in the sense of probability. 

The change by a monomer-event is as follows. If = Bfi = 0, that is, all nearest 
neighbors of site n are vacant, A adsorbs at site n, resulting in = 1. If Bft = 0 
but Afi 7 ^ 0, the adsorption attempt is neglected and no configuration change happens, 
which models strong repulsion between A’s. If Bfi 7 ^ 0 , we choose one site among 
S-occupied nearest neighbors of n with equal probability. Assume that site f is chosen. 
Then the adsorption-attempting A and B at site r react and form a ‘molecule’ AB 
which desorbs in no time. This event amounts to the change of af from 0 to —1. In 
any case including dimer-events below, a ‘molecule’ AB formed by reaction is always 
removed from the system immediately. 

A change by a dimer-event is as follows. If site m is not vacant, the adsorption 
attempt is neglected and the configuration remains the same. Assume that site m is 
also vacant. There are three cases which allow for a change of the conhguration. 

Case I. Afi Affi Bfi Bfjf 0. 

Case II: 7 ^ 0 and A^ = Bfh = 0; or Aff^ 7 ^ 0 and Afi = Bfi = 0. 

Case III: Aa 7 ^ 0, A^ 7 ^ 0. 

In all three cases, a dimer BB is hrst dissociated into two monomers which attempt 
to adsorb at site n and m, individually. In the case I, two BA adsorb at sites n and 
m, which gives CTff = arn = 0. To explain what will happen in the case II, we let s [f^ 
be the site with nonzero [zero] number of A’s in its nearest neighbors. Then, one B 
adsorbs at site r and another B reacts with one of A’s occupying nearest neighbors of 
s. This event results in cx^ = 0 and (7^= —1, where I is a randomly chosen site among 

§ Actually, a configuration in which all sites are occupied by particles does not allow any further 
dynamics. However, it is not considered an absorbing state, for no probability current to this 
configuration from any other state is present. 
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A-occupied nearest neighbors of s. In the case III, two B's react individually with 
randomly chosen A’s. That is, the conhguration becomes = — 1, where 

{rfin) is a randomly chosen site among A-occupied nearest neighbors of n (m). Except 
these cases, no conhguration change happens. 

After the adsorption attempt, regardless of whether it is successful or not, time 
increases by l/N^. 

For a given conhguration, we introduce a random variable A4 to be called ‘staggered 
magnetization’ (SM) as 

x = (1) 

i,j 

After many realizations with the same initial condition, we calculate the mean SM, 
m(t), and the (time-dependent) Binder cumulant, U(t), dehned as 

= = ( 2 ) 

where (.. .)t stands for the average at time t. At the critical point of the SBODT, m(t) 
decays to zero as where fd is the critical exponent for the order parameter, u 

is the correlation length exponent, and 2 ; is the dynamic exponent. We will estimate 
f3/{uz) from simulations and compare it with that of the Ising model. 

3. Simulation Results 

3.1. Preliminary : dynamic exponent of the two dimensional Ising model 

Since the main purpose of this paper is to hgure out whether the SBODT of the IMD is 
described by the critical exponents of the Ising model, we hrst present simulation results 
for the dynamic exponent of the Ising model at criticality, to be self-contained. This 
subsection may be read independently of other sections. In this subsection, a and m{t) 
should be understood as an Ising spin and the mean magnetization of the Ising model, 
respectively, and these should not be confused with the same notation for the IMD. 
The Ising Hamiltonian is 

L 

H J ^ ^ ^i,j T ^i,j+l) ) (3) 

where (Xjj is the Ising spin at site {i,j) and periodic boundary conditions are assumed. 
We remind that the critical point and energy per site at the critical point are exactly 
known as Kc = J/ksTc = ln(l -|- \/2)/2 and Ed J = — In 2/2 [I9]. For convenience, we 
set J = 1 and energy is measured in unit of J. 

We simulated the dynamics of the Ising model at the critical point, using single spin 
flip dynamics with the Metropolis algorithm. As an initial condition, a fully ordered 
state is used, that is, we set (Xjj = 1 for all at f = 0. We measure magnetization 
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m(t), fluctuation of magnetization V{t), and the energy per site E{t), defined as 



At the critical point, the asymptotic behaviors of these quantities are (see, for example, 

123 ) 


m[ 


i{t) = [1 + Bmt-^- + o{t-^-)\ , (7) 

e{t) = E{t) -Ec = [l + , (8) 

V{t) = [1 + , (9) 

where d is the dimensions of the system (d = 2 in this paper); /d = | and u = 1 are 
exactly known critical exponents (see, for example, |2T]); z is the dynamic exponent to be 
determined; A’s and S’s are constants; and y’s are exponents of the leading correction- 
to-scaling behavior which will be called leading corrections-to-scaling exponents (LCEs). 
To hnd the dynamic exponent, we investigate the effective exponent functions (EEFs) 
(6>1) 


^ \n[m(t)/m(t/b)] 1 ^ — 1 ^ / v ^ 

Sm{t] b) = ^ + o(r^™) 

In b QZ In 0 

^ ^ _ 1 _ + „(r-). 

\n[V{t)/V{t/b)] 7 ^b^^-1 ^ 

£.(*; b) = ' ^ Ld' " = -77--B.-T-i-*'*’ + »(«■*")■ 


In 6 


4z 


In b 


( 10 ) 

( 11 ) 

( 12 ) 


Since the LCE governs how the EEF approaches the asymptotic value, it is also 
important to hnd the value of the LCE. To this end, we use the method introduced 
recently [22l[23]. We numerically calculate the corrections-to-scaling functions (CTSFs) 
Qm{t-,b), Qe{t-,b), Qvit] b), dehned as 

m{f)m{t/b‘^) 


Qmit]b) = In 

Qe{t;b) = In 
Qvit-,b) = In 


m{t/by 

e{t)e{t/b'^) 

e{t/by 

V{t)V{t/b^) 


= Bm{bl,-l)h-^- + o{t-^-), 
= B,{b^^ - + o{t-^^), 


(13) 

(14) 

(15) 


V{t/bf 

Note that the knowledge of critical indices is not necessary to hnd y from Q. For 
convenience, we just drop indices like Q{t;b), £{t]b), x when we have to write EEFs, 
CTSFs, or LCEs collectively in the following. 

To estimate the critical exponents, we use the following strategy. First note that 
the correct value of y makes Q{t]b)/{b^ — 1)^ for any b he on a single curve Bt~^ in 
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Figure 1. Double-logarithmic plots of Q(t; b)/{b^ — 1)^ vs t for b = 10°'^ ~ 5, 6 = 10, 
and b = 10^'^ ~ 20 (left to right) with x = 1- Top three curves are for Qe{t\b) and 
bottom three curves for Qm{t; b). Asymptotically, curves for different 5’s lie on a single 
curve, suggesting Xm = Xe = 1- For comparison, we also plot B/x with B obtained 
from the htting of the effective exponents; see figure [2] and figure [S] 


the long time regime. Exploiting this feature, we adjust the value of x until double- 
logarithmic plots of Q{t]b)/{b^ — 1)^ as a function of t for different fe’s he on a single 
straight line. After Ending x, we plot S{t-,b) as a function of r = {b^ — l)f“^/ln6 for 
various 6’s. These curves should he on a straight line in the asymptotic regime, once 
X is correctly estimated. By extrapolating the straight line behavior of S{t]b), we can 
estimate the critical exponent. 

Now we present the simulation results. In simulations, the system size is L = 2^^, 
the maximum observation time ist = 10^, and the number of independent runs is 2 x 10^. 
We begin with the analysis of Qm and Qe- Figure [H depicts Q{t; b) /(6^ — 1)^ as functions 
of t for a few 6’s with x = 1- Since the asymptotic behavior does not depend on b with 
the choice of x = 1, the LCEs are estimated as Xm = Xe = 1- Note that this estimate 
is consistent with the fitting result by Nam et al |20j . 

We now analyze the EEFs. Figure [2] shows how Sm{t] b) behaves for sufficiently 
large t, when it is depicted against r = {b — l)/(tln6). Since EEFs are almost on the 
same line irrespective of b for r < 0.02, we conclude that the critical relaxation dynamics 
is well described by terms up to the leading correction for r < 0.02. Fitting of Emit] 10) 
for the region r < 0.015 using a linear function, l/(8z) — BmT with l/(8x) and to 
be fitting parameters, we get 

^ = 0.057 650(12), Bm ~ 0.114(2), (16) 

0.2 

or equivalently 

z = 2.168 26(45), 


(17) 
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(xlO-^) 



Figure 2. Plots of £m{t]b) +0.0575 vs t = [b— l)/(tln 6 ) for h = 10°-^, 10, and 10^-^. 
Note that the scale on the y axis is multiplied by 10“'*. In the region t < 0.02, the 
effective exponents lie on a single straight line. The fitting for r < 0.015 of the effective 
exponents gives 1 /( 80 ) = 0.057 650(12) and Bm = 0.114(2). Inset: Close-up view on 
the region r < 0.005 for 6 = 10 with the fitted linear function (straight line). 



T 


Figure 3. Plots of £e{t]h) vs r = ( 6 — I)/(tln 6 ) for b = 10°'^, 10, and 10*'^. The 
straight line is the result of the fitting £e{t\ 10) for the region r < 0.015, using a linear 
fitting function. Inset: Plots of as a function of 1/t for the magnetization (z^) 
and the energy density (Ze)- 


where the numbers in parentheses indicate the uncertainty of the last digits. We use 
this convention throughout the paper. 

Figure |3] depicts the behavior of Se{t] b) for some 6’s. Similar to Sm, all curves 
for r < 0.02 lie on a single straight line. Due to the statistical noise of the data, the 
estimate of z is less accurate than (1T7|) . Still, we get a consistent result within error 
























Ising criticality of the IMD model 



Figure 4. Log-log plot of L?V{t) vs t. A line segment with slope 7({Az) ~ 0.8071 is 
for a guide to the eyes. Inset: Plot of £v vs 1/t for b = 10. 


bars. For comparison, we define effective exponents for z as 

1 ... 1 


Zmit) = 


Ze{t) = 


(18) 


8SUt]ioy ^e(t;io)’ 

which are collectively called The inset of figure [3] compares these two effective 
exponents, which shows an agreement of the limiting values of two curves within errors. 

Figure m depicts L‘^V{t) as a function of t. Since V(t) is the fluctuation of m(t), 
it is not surprising that V{t) is much noisier than m{t). The noisy behavior of V{t) 
becomes conspicuous when is depicted (see the inset of figure 0]). Since is quite 
noisy compared to Sm and Sg, the error of the estimate of the leading scaling exponent, 
7/(4^), solely from is larger than before. From we estimate 7/(4^) = 0.81(1) 
which is, of course, consistent with the estimate flT7)) within errors. 

As we have seen, statistical noise is minimal when the relaxation of magnetization 
is investigated. Thus, we conclude (d/ivz) = 0.057 650(12) and, accordingly, 2 ; = 
2.168 26(45) for the single spin flip dynamics of the two dimensional Ising model with 
the Metropolis algorithm. 


3.2. Interacting monomer-dimer model 

In this section, we present the simulation results for the IMD. We begin with analyzing 
the EEF £m{t) and the CTSF Qm{t) for the SM, defined similarly as ffTOjl and flT^ . 
respectively. Since we do not know the critical point a priori, £m{t) will be used to 
find the critical point, by exploiting the fact that £m(t) should veer up (down) as t gets 
larger if the system is in the ordered (disordered) phase. 

The initial condition of simulation is that all odd sites are vacant and each even 
site can be either A-occupied with probability thq or vacant with probability 1 — tuq. 
With this initial preparation, the SM at time t = 0 is m(0) = mo. As in [16], we set 
mo = 0.7. Figure [5] depicts £rn with 6 = 10 as a function of for p = 0.638 59 
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Figure 5. Plots of the effective exponents Em of the SM against for p = 0.638 59 
and 0.6386. Here, we are using b = 10. Inset: Log-log plots of Qml{b^ — 1)^ vs t with 
X = 0.75 for b = 10°'^, 10, 10^'^ (left to right). The straight line with the slope 0.75 is 
also drawn for comparison. 


and 0.6386. Here, the system size is L = 2^^, the number of independent runs is 1400, 
and the system evolves up to t = 10®. The effective exponent for p = 0.638 59 [0.6386] 
approaches the Ising value (1T6|1 then veers down [up] for large t, suggesting that the 
critical point is Pc = 0.638 595(5). This observation shows that the critical exponent of 
the IMD is consistent with the Ising value, contrary to the claim in |16] . 

The inset of hgure |5] shows the behavior of Qm/{b^ — 1)^ at p = 0.638 59 with 
X = 0.75 for some values of 6 on a double-logarithmic scale. Since Qm is measured at 
the disordered phase, the curves should eventually veer down [23] as can be seen at the 
tail of the curves in the inset. Still, the power-law region is observable and we conclude 
that the CSE is about 0.75. Notice that the LCSE of the IMD is smaller than that of 
the Ising model (see hgure [T|). Thus, to hnd the correct value of critical exponents for 
the IMD requires longer evolution time than the Ising model, which is also observed 
in [6] for the 2DIM. 

The study of £m suggests that the order-disorder transition of the IMD is indeed 
described by the critical exponents of the Ising model. To support our claim further, we 
also study the behavior of the Binder cumulant U{t). As in [16], the initial condition 
for the study of the Binder cumulant is the fully vacant state without any A and B. 
In hgure |6l we showed the behavior of U(t) at p = 0.638 59 (top red curves) and 
p = 0.638 55 (bottom blue curves) for L = 2® (left two curves) and 2^® (right two 
curves). The number of independent runs for L = 2® and p = 0.638 55 (p = 0.638 59) 
is 40 000 (60 000) and the number of independent runs for L = 2^® and p = 0.638 55 
(p = 0.638 59) is 2000 (3000). At p = 0.638 55 which is the estimated critical point in 
[T6] . we also observed that U(t) approaches around 0.6 for L = 2® as in [16]. However, 
U{t) for L = 2^® approaches 0.58, indicating that the system with p = 0.638 55 is in the 
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Figure 6. Plots of the time dependent Binder cumulant against t on a semi-logarithmic 
scales for L = 2® (left two curves) and L = (right two curves) at p = 0.638 55 (top 
red curves) and p = 0.638 59 (bottom blue curves). For comparison, a straight line 
indicating the critical Binder cumulant of the Ising model (0.611) is also drawn. 


disordered phase. This observation is consistent with onr estimate of the critical point. 

As in the 2DIM stndied in [6], analyzing the Binder cnmnlant is not an efficient 
method to hnd the critical point of the SBODT for models with two symmetric absorbing 
states. Even thongh the system with p = 0.638 59 is in the disordered phase, the Binder 
cnmnlant increases up to L = which might lead to a wrong conclusion. Since the 
Binder cumulant at p = 0.638 59 is quite close to the Ising critical value for L = 2^°, 
to hnd the correct critical point using Binder cumulant requires the system size to be 
larger than L = 2^^. 

4. Summary and Conclusion 

We numerically studied the two-dimensional interacting monomer-dimer model (IMD), 
focusing on the symmetry breaking order-disorder transition. Relaxation dynamics of 
the ‘staggered magnetization’ around the critical point was analyzed. The analysis 
of the corrections to scaling function showed that the two-dimensional IMD model 
has stronger corrections to scaling than the Ising model. We found that the critical 
point of the IMD is Pc = 0.638 595(5) and the relaxation dynamics at the critical 
point is consistent with the critical relaxation of the Ising model which is estimated 
as (d/iyz) = 0.057 650(12) in section [XU Thus, we concluded that the order-disorder 
transition in the two dimensional IMD shares criticality with the Ising model. As a 
hnal remark, we would like to emphasize that due to strong corrections to scaling, the 
estimate of the critical point using the Binder cumulant becomes unreliable unless the 
system size is larger than 2^*^. 
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